# disable printing of GC info
kernelopts(printbytes=false);

# modify here for different parity-check matrices.
#read "H1008.mpl":
#read "H19982.mpl":
read "H4376.mpl":

# prune the matrix to make it linearly independent
# FIXME: need to write a procedure to automate this process
H := LinearAlgebra:-DeleteRow(H, 1):

# generate parity-check matrix with dimension M * N
(M, N) := LinearAlgebra:-Dimension(H);
K := N - M;

read "util.mpl";
read "encoder.mpl":
read "scrub.mpl":

t := columnWeight(H);

# generate matrix with dimension K * N
G := computeG(H);

# rate >= R, if H is linearly independent, then rate = R.
R := evalf(K/N);

# maximum iteration number
maxiter := 32;

# seed the randomizer with system clock.
Seed := randomize():

# generate a random sequence of bits of length K as the transmitted
# codewords.
genRandomMessage := proc()
global G, K;
    encodeLDPC([seq((rand() mod 2), i=1..K)], G);
end proc:

# verify the correctness of the generator matrix G computed above
verifyG := proc()
global H;
    evalb(add(i, i=LinearAlgebra[Modular]:-Multiply(2, H, Vector(genRandomMessage()))) = 0)
end proc:

# cell level
cellLevel := 4;

# cells per segment
# 9 for 1998 length block
# there are 9 segments in total
#cellsPerSegment := 111; # for H19982.mpl, i.e., 9 segments
cellsPerSegment := 547; # for H4376.mpl, i.e., 4 segments

# time to stop measurement
timeToStop := 6 * 3 * (10**6);

# scrubbing gap for periodical scheme
scrubGap := 5 * (10**3);

# scrubbing gap for basic scheme
#scrubGapBasic := scrubGap * iquo(iquo(N, log[2](cellLevel)), cellsPerSegment);
scrubGapBasic := 6 * (10**4);

# performance counters for recording LDPC decoding performance in periodical scrubbing.
averageBitErrorRates := [seq([0, 0], i = 1..(timeToStop/scrubGap))]:
averageIterations := 0:
averageSuccessfulRate := 0:
averageLastBER := 0:
averageLastSuccessfulRate := 0:

# directory and filename which stores the experimental data.
savedir := "/Users/yli/Desktop/Experiments/data/";
filename := cat(savedir, cat( cat(StringTools:-FormatTime("%m-%d-%H-%M-%S"), ".mpl")));

# update and save the values of global performance counters to file.
updateAndSaveData := proc(stats, trials, succno)
global averageBitErrorRates,averageIterations, averageSuccessfulRate, averageLastBER, averageLastSuccessfulRate;
global filename;
    averageBitErrorRates := averageBitErrorRates/trials;
    averageIterations := stats[1];
    averageSuccessfulRate := stats[3]/stats[2];
    averageLastBER := stats[6];
    averageLastSuccessfulRate := evalf(succno/trials);
    save averageBitErrorRates,averageIterations, averageSuccessfulRate, averageLastBER, averageLastSuccessfulRate, filename;
end proc:

# simulate scrubbing, returning the number of faulty bits.
# blkno: number of trials
# opt: optimization swtich to determine if we use error log before
# ECC.
# choice: scrubbing scheme selection:
# 1 for our staged scrubbing
# 2 for the basic periodical scrubbing
simulateScrub := proc(blkno, opt, choice)
local i, stats, stat, succno;
global cellLevel, cellsPerSegment, timeToStop, scrubGap, scrubGapBasic, maxiter, H, N;
global averageBitErrorRates;

    succno := 0;
    stats := [0, 0, 0, 0, 0, 0];
    if choice = 1 then
        print("Using staged scrubbing");
    elif choice = 2 then
        print("Using basic scrubbing");
    else
        error "unrecognized choice of scrubbing scheme";
    end if;
    for i from 1 to blkno do
        myprint("Number of trial: ", i);
        if choice = 1 then
            stat := periodicalScrub(genRandomMessage(), cellLevel, cellsPerSegment, scrubGap, timeToStop, maxiter, H, opt);
        else
            stat := basicScrub(genRandomMessage(), cellLevel, scrubGapBasic, timeToStop, maxiter, H);
        end if;
        stats := stats + [seq(stat[i], i=1..6)];
        averageBitErrorRates := averageBitErrorRates + stat[7];
        if stat[7] = 0 then
            succno := succno + 1;
        end if;
    od;
    stats[1] := evalf(stats[1]/blkno);
    stats[6] := evalf(stats[6]/(blkno * N));
    updateAndSaveData(stats, blkno, succno);

    print("Final statistics of the advanced scrubbing simulation:");
    # 1. average termination iteration
    # 2. no of segments which are scrubbed
    # 3. correct block no
    # 4. detected block no
    # 5. undetected block no
    # 6. bit error rate.
    print(stats);
    stats
end proc:

# simulate naive scrubbing using LDPC code
simulateNoScrub := proc(blkno)
global cellLevel, timeToStop;
global H, M, N, K, t, R, G, maxiter;
local i, trans, recv, corrected, decodedResult, successflag;
local errbits, totalErrBits, totalErrorBlocks, totalIters;
local correctblkno, undetblkno, detblkno, stat;
    # initializing counters
    correctblkno := 0;
    undetblkno := 0;
    detblkno := 0;
    totalErrBits := 0;
    totalErrorBlocks := 0;
    totalIters := 0;

    # start decoding
    for i from 1 to blkno do
        myprint("Now sending block:", i);
        trans := genRandomMessage();
        # a pair is returned (flag, decodedMessage)
        # currently using sum-product algorithm (not in log-domain)
        decodedResult := noScrub(trans, cellLevel, timeToStop, maxiter, H);
        successflag := decodedResult[1];
        totalIters := totalIters + decodedResult[2];
        corrected := decodedResult[3];
        errbits := countErrorBits(trans, corrected);
        if errbits = 0 then
            if successflag then
                correctblkno := correctblkno + 1;
            else # this branch is not quite possible to happen unless a bug.
                error "warning, a bug happens";
            end if;
        else #if errbits > 0 then
            totalErrorBlocks := totalErrorBlocks + 1;
            totalErrBits := totalErrBits + errbits;
            if successflag then
                # then this is an undetected error
                undetblkno := undetblkno + 1;
            else
                detblkno := detblkno + 1;
            end if;
        end if;
    od;

    # compute statistics
    # bit error rates, block error rates, average iteration number
    print("Simulation statistics:");
    stat := [evalf(totalErrBits/(N*blkno)), evalf(totalErrorBlocks/blkno),  evalf(totalIters/blkno), correctblkno, undetblkno, detblkno];
    print(stat);
    stat
end proc:

#time(simulateNaiveScrub(100));
#time(simulateScrub(100, true));
